function nidat = fsb_filter_spatial_volume(indat,filtdat,Tempvar)

% FSB : Spatially filter image data scanwise
%
% EXAMPLE:
% idat = fsb_filter_spatial_volume(idat,sandbox)
%
% INPUT:
% idat : 4D-matrix containing image data
% filtdat: struct containing filter information
% Tempvar: struct containing ROI information
%
% OUTPUT:
% idat: spatially filtered image data
% 
% CALLED BY:
% FSB
%
% NOTES:
% Filter used is a Butterworth filter
% Fc : Filter cutoff in Hz
% Fs : Sampling frequency in Hz
% x  : Data, vector or matrix (works on columns then)
% type: can be 'high' or 'low' for a high- or lowpass filter or 'stop' for
% a bandstop filter
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check input
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<2;
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Load default values
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    filtdat.filt_freq = [5 100];
    filtdat.filttype = 'stop';
    filtdat.roi =5;
end

try
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Look for filter type input
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    f_type = filtdat.filttype;
catch
    disp('Filter type not defined, reverting to Bandpass filter with standard values')
    f_type = 'Stop';
    filtdat.filt_freq = [1/filtdat.filt_freq_high 1/filtdat.filt_freq_low];
    filtdat.filttype = 'Stop';
end

f_name = f_type;

if strcmp(f_name,'Stop')==1;
    filtstring = ['Bandpass filtering, High cutoff: ' num2str(filtdat.filt_freq(1)) ' s , Low : ' num2str(filtdat.filt_freq(2)) ' s'];
else
    filtstring = ['Filtering with ' f_name 'pass filter, threshold : ' num2str(filtdat.filt_freq(1)) ' s'];
end

dim = size(indat);
if size(dim,2)<4;
    dim(4) = 1;
end
nidat = indat;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Apply filter
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if strcmp(f_type,'Spectrum')~=1;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Set up Filter
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if size(filtdat.filt_freq) ==1
        Fc = 1/filtdat.filt_freq(1);
        Fs = 1/filtdat.filt_freq(2);
        ord = 4;
        filt = Fc*2./Fs;
        disp(['Filter is : ' num2str(filt)]);
        if max(filt)>1 || min(filt) <0
            disp('Filter size not possible, aborting')
            return
        end
        [B,A] = butter(ord, filt,f_type);
    else
        Fc=[1 1]./filtdat.filt_freq;
        Fs = [1/2 1/2];
        ord = 2;
        filt = Fc*2./Fs;
        disp(['Filter is : ' num2str(filt)]);
        if max(filt)>1 || min(filt) <0
            disp('Filter size not possible, aborting')
            return
        end
        [B,A] = butter(ord, filt);
    end
    h = waitbar(0,filtstring);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Go over volumes
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    for scannum = 1:dim(4)

        idat3 = indat(:,:,:,scannum);
        mean_idat2 = mean(idat3(:));
        idat3 = idat3(:);
        det = filtfilt(B,A,idat3);
        idat3 = reshape(idat3,dim(1),dim(2),dim(3));

%         %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%         % 2-D version, in testing
%         %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 
%         idat3 = shiftdim(idat,3);
%         waitbar(5/10);
%         idat3 = reshape(idat3,size(idat3,1),[]);
%         waitbar(6/10);

%        det = filtfilt(B,A,idat3);

%         waitbar(8/10);
%         idat = reshape(det,[dim(4),dim(1),dim(2),dim(3)]);
%         waitbar(9/10);
%         idat = shiftdim(idat,1);

        idat3 = idat3+mean_idat2;
        mean_idat = mean(idat3(:));
        disp(['Mean idat value after = ' num2str(mean_idat)]);
        nidat(:,:,:,scannum) = idat3;
    end

    close (h);
    disp('Filtering done')

else

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Calculate raw data average and their frequency spectrum
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    Fs = 1/filtdat.TR;
    % Set parameter for mean volume calculation to zero to speed up things for
    % testing
    allvol = 1;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % calculate ROI data
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    try
        idat2 =nidat(Tempvar.slice_n(1)-Tempvar.roi:Tempvar.slice_n(1)+Tempvar.roi,...
            Tempvar.slice_n(2)-Tempvar.roi:Tempvar.slice_n(2)+Tempvar.roi,...
            Tempvar.slice_n(3)-Tempvar.roi:Tempvar.slice_n(3)+Tempvar.roi,:);
    catch
        disp('ROI size exceeding volume, aborting')
        return
    end
    idat2 = shiftdim(idat2,3);
    idat2 = reshape(idat2,dim(4),[]);
    idat2 = single(idat2);
    iROI = mean(idat2,2);
    idat2 = idat2-repmat(mean(iROI(:)),size(idat2));

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Calculate parameters for image display
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    L = length(iROI);
    NFFT = 2^nextpow2(L); % Next power of 2 from length of L
    f = Fs/2*linspace(0,1,NFFT/2);

    h = waitbar(0,'Analyzing spectrum for ROI...');

    Y2 = abs(fft(idat2,NFFT)/L);

    ROI_spect = mean(Y2,2);
    close (h);

    if allvol ==1;

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Calculate center of volume data
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        idat3 = nidat(round(dim(1)/4):round(dim(1)/4*3),...
            round(dim(2)/4):round(dim(2)/4*3),...
            round(dim(3)/4):round(dim(3)/4*3),:);
        idat3 = shiftdim(idat3,3);
        idat3 = reshape(idat3,dim(4),[]);
        idat3 = single(idat3);
        imean = mean(idat3,2);
        idat3 = idat3-repmat(mean(imean(:)),size(idat3));

        h = waitbar(0,'Analyzing spectrum for central part of image...');
        Ya = abs(fft(idat3,NFFT)/L);
        mean_spect = mean(Ya,2);
        close (h);
    end

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Plot raw data average and their frequency spectrum
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if allvol ==1;
        figure(601); cla; set(gcf,'Name','Volume center signal time course');

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Plot mean volume intensity time course
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        subplot(2,1,1);
        plot(imean);
        title('Mean signal time course')
        xlabel('time(volumes)')
        ylabel('Signal intensity')

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Plot single-sided amplitude spectrum.
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        subplot(2,1,2)
        plot(f,2*mean_spect(1:NFFT/2));
        title('Single-Sided Amplitude Spectrum of mean volume intensity(t)')
        xlabel('Frequency (Hz)')
        ylabel('|Y(f)|')

    end

    figure(603); cla; set(gcf,'Name','ROI signal time course');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Plot mean ROI intensity time course
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    subplot(2,1,1);
    plot(iROI);
    title('Mean ROI time course')
    xlabel('time(volumes)')
    ylabel('Signal intensity')

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Plot single-sided amplitude spectrum.
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    subplot(2,1,2)
    plot(f,2*ROI_spect(1:NFFT/2));
    title('Single-Sided Amplitude Spectrum of mean ROI intensity(t)')
    xlabel('Frequency (Hz)')
    ylabel('|Y(f)|')

end
end




